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Abstract 

We show that the boundary curves (profiles) in R 2 of the generalized projections of a body in R 3 
uniquely determine a large class of shapes, and that sparse profile data, combined with projection 
volume (brightness) data, can be used to reconstruct the shape and the spin state of a body. We also 
present an optimal strategy for the relative weighting of the data modes in the inverse problem, and 
derive the maximum compatibility estimate (MCE) that corresponds to the maximum likelihood or 
maximum a posteriori estimates in the case of a single data mode. MCE is not explicitly dependent 
on the noise levels, scale factors or numbers of data points of the complementary data modes, and 
can be determined without the mode weight parameters. We present a solution method well suitable 
for adaptive optics images in particular, and discuss various choices of regularization functions. 
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1 Introduction 

Constructing the shape model of a three-dimensional object or a surface is often based on images ob- 
tained at various viewing geometries. This is the standard case in human and robot vision as well as 
in cartography. When individual points on the surface can be identified in different images, one can 
directly solve for their position vectors and use this as a top-down basis for stereographic mapping or 
model construction. In many cases, however, the model construction reduces to an inverse problem rather 
than a cartographic one as the image information is based purely on the projections of the target in the 
viewing directions. When surface illumination and other effects are taken into account, we talk about 
generalized projections |14j to distinguish these from simple shadow-like projections or silhouettes. As 
discussed in [14j . there are various types of such projections, ranging from the volume- like quantity of 
integrated brightness Lei (generalization of the area of a shadow on a projection screen) to resolved 
images I £ I 2 x K (for one wavelength) . 

In this paper, we consider the case where images Z(w,wo) obtained at various viewing and illumination 
directions w, ojq G S 2 are available, but the infomation in these images is only contained in the boundary 
curves between the dark background or a shadow and the illuminated portion of the target surface. This 
situation is typical for faraway objects in space for which low-resolution images are available via large 
telescopes through adaptive optics (AO) [19] or other deconvolution and image processing techniques. 
Due to the deconvolution process, the actual brightness levels of the pixels in these images tend to portray 
artificial and exaggerated features, so they are usually less reliable than profile contours, i.e., the locations 
of the light/dark boundary pixels [U (3J. 



As the coverage of viewing geometries is seldom wide enough to enable a full reconstruction of the model 
from images alone, we also consider the possibility of augmenting the image dataset with a set of measured 
brightnesses L(uj,ujq) of the target at various observing geometries. Subsets of these, measured within 
some time intervals, are called lightcurves. As discussed in jTU]-[Il] (see also references therein), a global, 
usually convex, model of the target can be obtained from a large enough set of L{u),ujq). The possibility 
to use images I both serves to reconstruct more details in the model and to use a combined dataset for 
successful modelling when neither the available L nor T are sufficient alone. 

The ratio of profile contour pixels to all pixels of the target is approximately 4/D, where D is a typical 
diameter of the target in pixels. For the largest few asteroid AO targets, D is around 30. When D is 
less than about 10, the location accuracy of the border pixels is not necessarily very much better than 
the brightness accuracy of the pixels, but on the other hand there is no particular loss in the number 
of information points when only borders are used. For higher D, border points are a smaller subset of 
available pixels, but now their location accuracy is far better than the brightness accuracy of all pixels 



The profile contours of the AO images are obtained as a solution of a separate inverse (or imaging) prob- 
lem, where an approximation of the atmospheric point-spread function (PSF) is first used to deconvolve 
the raw image (with, e.g., connectedness of the processed image as a constraint), and the contours can 
be separately modelled with wavelet techniques [U [3J Q35] . This independence from the actual model of 
the target is advantageous in the sense that the assumptions and inevitable deficiencies of the model 
(particularly in the adopted light-scattering model on the surface [15]) do not affect the outcome of the 
AO image processing. On the other hand, from the methodological point of view, all information should 
usually be employed simultaneously when solving an inverse problem, so another approach would be to 
use the raw AO image data and the approximated PSF directly in model construction without separate 
image deconvolution. However, in practice it appears that the profile curve extraction procedure, in par- 
ticular, retains valuable independent information [2,3, and the model deficiencies affect the fit deviation 
between the predicted and observed model profile curves less than they affect that between the full model 
and AO images. What is more, below we show that the profile curves convey almost as much information 
on the shape as the full images, so we can conclude that the two-step inversion of AO (and brightness) 
data is well justified. 

The paper is organized in sections in the following manner. In Section 2 we study the information content 
of profile contours and show uniqueness results for the inverse problem of reconstructing shapes from 
these. Section 3 deals with the posing of the inverse problem and the choice of regularization functions. 
In Section 4, we discuss the weighting and maximum compatibility estimates for inverse problems with 
multiple data modes, and examples of the use of real brightness data L and images I are presented in 
Section 5. Section 6 sums up. 

2 Generalized profiles: uniqueness results 

In this section, we define the concept of generalized profiles as boundary curves corresponding to gen- 
eralized projections, and show that a large class of shapes is reconstructable from these. While many 
practical procedures for shape-from-profiles (also known as volume carving) and shape-from-shading are 
well known in, e.g., computer vision (see |21j and references therein) and cartography (clinometry) , some 
of their geometric characteristics and the properties of the corresponding shape classes and inverse prob- 
lems discussed here have not been previously stated in the mathematical literature, to the best of our 
knowledge. 

Let us first look at classical profiles defined by one direction uj G S 2 , i.e., loq — u>. The projection 
V(oj,B) G K 2 of a compact set 8gl 3 (a set of points on closed surfaces) in the direction uj £ S 2 maps 
x e B ->• x e V: 
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where < -d < tt,0 < ip < 2tt are the polar coordinate angles defining 

uj = (sin^cos tjj, sin#sin^>, cosi?), 

and Ri(a) is the rotation matrix corresponding to the rotation of the coordinate frame through angle a 
in the positive direction about the i-axis. Thus, e.g., R z (4>) 1S 

cos <p sin 

Rz(4>) = ( -sin0 cos<-6 ) . (2) 
1 

For "3 = 0, Xi = Xi, i = 1, 2, and for # = n, x\ = —xi,X2 = xi. Only a half of S 2 is needed for defining 

uj as 

Xi{ui) = (-l) l Xi(-u>). 

For some definitions below, we need to give the three-dimensional position x x (u}) S R 3 of the planar 
points h: 

.r^U) =R,( -i ) R„ ( - | X! | . (.T) 



x 2 

Definition 1. The profile Vd(w,B) of £> in the direction uj is the boundary of its projection V(cj,B): 

r d (u,B) = &P(LJ,B) 

(the specific notation Vq is used for emphasis). More specifically, Vo^^B) is the mapping x — > x from 
those points x £ B for which, for all lines x + sui parametrized by s, 

V a {uj, B) = {^|Vs : 36, Bp > : ¥(x + suj; 0, g) £ B, < g < p} , (4) 

where 0, g) denotes a parallel transport, perpendicular to the line £, of C by the amount g in the 
direction 6* e 5* 1 in some system around C. 

Definition 2. The cylinder continuation (CC) C(uj 7 S) of a set 5 of points x £ M. 3 is the set of points in 
M 3 given by 

C(u>,S) = ja; + su> x € «S; — oo < s < oo|. (5) 

Definition 3. The profile hull H e R 3 is the bounding surface of the set of points formed by the 
intersection of the CCs of the projections V in some directions Wj, i = 1, ...,iV, corresponding to 
measured profiles VdiiOi) € R 2 : 

«({Pe(wi)|£i}) = 0n c ( w *' 5< )' -Si = {Mwi)|xe7'(a;i)}. (6) 



Remark. For practical purposes and clarity, we assume H to be constructed such that it has a closed 
surface as a boundary, though this is not strictly necessary in its definition (we could define % as set of 
points rather than its bounding surface). Thus B and Ji are in the same object class. Similarly, we use 
the concepts of projection and profile (or body and surface) somewhat interchangeably when the meaning 
is obvious. 

Many bodies can be reconstructed to arbitrary accuracy with profile hulls. A convex body is already 
determined by its Ti. constructed with a full coverage of the directions ui confined to any plane in R 3 . In 
real profile measurements, however, the position of the profile in the xr-plane is usually arbitrary, i.e., the 
profile is determined up to a translation xo of the profile plane origin. Then data restricted to planar 
directions are not necessarily sufficient even for convex bodies: the profile offsets xq^ of a convex body 
are not always uniquely defined by the profiles Vd(u>i) via the profile hull H when u>i are confined to a 
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plane. This is simple to illustrate by considering curves in R 2 and their projections in R at directions in 
S 1 . 

If the curvature function C : S 1 — > R of a closed convex curve in R 2 is given by a real- valued non-negative 
Fourier series 

C(tp) = V c n e inv > 0, n > 0, 



where tp denotes the direction of the outward normal of the curve, the projected width w((p) of the contour 
in that direction is 

fTr/2 



where 



and for n/1 



w(<p) = / C(^ + v)cos^# = 3?V Cne in<p I n 

-tt/2 



tt/2 

cos n - cos "0 cfaA 

•rr/2 



n — 3, 5, 7, . . . 
2/(n 2 -l) n = 2,6,10,... 
2/(1 -n 2 ) n = 0,4,8,... 

and ci = since 7i ^ and we must have w(ip) = w((p + 7r). Thus w(ip) carries no information on the 
odd- valued n-coefficients of the curvature function C (<p) that uniquely defines the shape of the curve (cf. 
[TU] for convex surfaces in R 3 ). Thus, for cylindrical convex surfaces in R 3 , the observed profiles in the 
symmetry plane can be made to correspond to any odd parts of C((p) with suitable chosen offsets xro- A 
typical case is that of shapes mimicking a circular cylinder with constant C in the symmetry plane: for 
example, if we change kq = to 

4%) = i?(-ir div i [A C os(^mod^ -I)- 1], (7) 

the reconstructed shape is a cylindrical Reuleaux triangle. This degeneracy occurs since for projections 
R 2 — > R the volume of a profile is equivalent to the profile itself up to an offset. An additional profile 
from a direction perpendicular to the plane resolves the degeneracy. 

Thus, in general, we need data at full ui € S 2 for a unique reconstruction of a body when the profile 
offsets are not known. The principle in the reconstruction via the profile hull % is the requirement that 
% must be consistent with the observed profiles, i.e., the profiles of the constructed % must be identical 
to the observed ones: 

V d [ij i ,H({V d (u; j )\? =1 })]=P d (u k ) 
(otherwise the volume of the intersection defining % is not maximal). 

Let us denote by the complete profile hull %c the profile hull for which uj covers all of S* 2 . The complete 
profile hull Hc(B) of a body B is the envelope of those of its tangents that do not intersect B anywhere. 
Then we can define a class of surfaces that includes convex ones but extends far into nonconvex surfaces: 

Definition 4. Tangent-covered bodies (TCBs) are bodies that are their own complete profile hulls: 
B = Wc(B). Thus, each surface point x G B of a TCB is mapped at least to one Pa(cj). 

TCBs include a large variety of nonconvex surfaces or sets of them: for example, a body consisting of 
two separate spheres is a TCB. While convex bodies C are reconstructable from the volumes of their 
generalized projections [TUJ [14] and are defined by having no tangents intersecting the body, TCBs T are 
reconstructable from profiles and are defined by there being at least one tangent at each surface point 
not intersecting the body elsewhere. 

Let us now generalize the concept of profile in the same way as projections. This leads to a shape class 
G, larger than TCBs T, that can be reconstructed from generalized profiles: 

CcTcG. 
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When we consider the directions (uj, ujq) in S 2 x S 2 , the region both visible and illuminated on B is given 
by [Hill] 

A + {uj,uj q ;B) = A + {uj;B)^A+{uj q ;B), (8) 

where 

A+(uj;B) = ja: G B (i/(a;),u;) > 0; Vs > : x + £ , (9) 

where f (a;) is the unit surface normal at x. The projection V of the boundary dA+ is now the generalized 
profile: 

Definition 5. The generalized profile of the body 23 in the direction uj and at illumination direction ujq 
is 

dP[u,A+(u,wo;B)]=V[w,dA+(u,uo]B)]. (10) 

The shape class Q is not as straightforward to define as C and T . We can, however, prove configurations 
allowing unique shape determination that illustrate its extension from T. 

Theorem 1. Assume that we know some parts K, of a body B from profile measurements. There exist 
configurations in which unknown parts U of B not determinable from profiles can be uniquely determined 
from the generalized profiles of B by using the shadow boundaries of K. on 11. 

Proof. Assume that He (B) is defined, and that it contains a planar section K,, and that B contains in 
this region an unknown concavity U (corresponding to He \ B). Also, assume that all points of U are 
seen from the viewing direction w_L/C, and that the illumination direction ujq lies in a plane -L/C, with 
9 = £(u)q, /C). Then the planar edge curve of the concavity U can be determined when 9 — > 0: 

dU = lim dS(9) := dS , 

0-s-O 

where dS denotes the projection of the shadow boundary in the direction u on the plane JC. At various 
< 9 < 7r/2, we can measure the shadow boundary projections dS(9), and thus extract the projection 
dS(9) of the shadow boundary inside U: 

dS{9) = dS{9) \ dS(9), dS(9) := dS(9) n dS - 

Then the envelope in K 3 of the intersection curves of the cylinder continuations of dS(9) in the directions 
of u> and loq 

C[u,d§(9)]nC[u o ,d§{0)] 
uniquely constructs the surface of the concavity U (when U is suitably regular). □ 

Theorem 2. There exist configurations in which unknown parts U of B can be uniquely determined by 
using their shadow boundaries on the known part K,. 

Proof. Let B be a combination of a TCB £ and any surface T> (£ n T> = 0) that can be determined using 
profiles in the directions uj for which the profile intersection 

Q(oj) :=V(oj,V)nV(u>,£) 

vanishes, Q — (the whole of T> is in the known part JC). The unknown parts are assumed to be on £ 
(they cannot be determined using the above uj). Now the parts of profiles of £ that merge with V(oj, V) at 
some uj, i.e., Q(u>) ^ 0, are represented as shadows on T> that we assume we can see from some directions 
uj' . The full or partial profiles of &P(uj,£) for which dV(uj,£) D Q = can be determined as usual, and, 
with a known D, the remaining parts dV{yj, £)flQ^0 can be determined from the shadows on T>. The 
intersection 

dW = C[oj',dS p (uj')}nV, 

where C[uj' , dS p (u)')] denotes the cylinder continuation corresponding to the observed projection of the 
shadow boundary of £ on T> in the direction uj' , can be used to determine the projection 'P(<9yV, uj), which 
completes the missing parts of the needed profiles. Now we have constructed the set of full profiles of £ 
at all uj £ S 2 , so £ can be constructed as it is a TCB. □ 
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Continuing in a similar manner, we can construct more complex variations of the two cases above to 
explore the shape class Q. In practice, directions (oj,luq) seldom cover S 2 x S 2 extensively or densely, 
so the shape is reconstructed within some resolution (discretization degree of the model) and a priori 
assumptions, as discussed below. 

3 Inverse problem 

Let us now consider the inverse problem of determining the shape and spin state of a body B from some 
measured generalized profiles dV[u>i, A+(uii, woi! B)] and the volumes L(ojoi, u>i) of generalized projections. 
We present a method that is suitable for typical ground-based astronomical data, i.e., the profiles are 
only obtained at restricted geometries and their resolution level is not high. When a dense coverage of 
geometries and high resolution are available (e.g., space probe missions), direct methods of computer 
vision and cartography are usually applicable. 

Our goal is to construct a total goodness-of-fft measure Xtot 

Xtot = xl + ^dxl + ^Rg(P), (ii) 

where L denotes lightcurves, d generalized profiles, and R regularizing functions g(P), where P £ W is the 
vector of model parameters. Regularization is discussed at the end of this section, and the determination 
of the weights A in section 4. We note here that the additional g(P) make xtot pseudo-x 2 as it no longer 
describes an underlying (assumed) Gaussian probability distribution (though g(P) may be x 2 -like in their 
functional form). When probability densities such as a posteriori distributions are constructed from Xtot; 
one can assume ^-distributions (of the form e~ cx ) only for the data components, and other suitable 
(prior) distributions for the regularization components |17) such that the maximum of the a posteriori 
distribution occurs at argminx 2 t(-P)- 

Throughout this paper, we do not include the conventional l/<5 2 -factor in ^-forms, where 8 is the 
expected (Gaussian) error variance (noise level), since 8 is seldom known exactly, and it does not affect 
the determination of our point estimates which is the goal of this paper. Suitable parameters for Gaussian 
or other distribution widths can be inserted separately whenever we want to construct a distribution. 

The volumes of generalized projections are also called total or disk-integrated brightnesses [TOl [Til [14] : 

L{uq,w) = I R{x;u>o,u>)(ui,i'(x)) da(x) = / i?[P _1 (oj, A+, h)\ cjo, ui]d 2 x, (12) 

where v{x) and da{x) are, respectively, the outward surface normal and surface patch of £>, R(x;uJa,w) 
describes the intensity of scattered light at the point x on the surface, P _1 (w,yl + , x) is the point in A+ 
corresponding to the projection point x, and d 2 x is the surface patch of the projection V . In its basic 
form, 

xl = £> (ob, W«0 - £ (mod W^)] 2 (13) 

i 

(assuming a constant noise level; see [11] for modifications and variations of this). L-data on S* 2 x S 2 
uniquely determine a convex body and the solution is stable [TU1 [H] , but L-data do not carry information 
on nonconvexities in most realistically available S 2 x S 2 geometries in practice [3] . Such information must 
be provided by AO or other techniques. 

For many typical AO targets, the generalized profiles are starlike due to the proximity of lo and wo and 
some regularity of the target shape at global scale [2J [3] [T5]. Then we can write x% by considering, for 
each profile i, their observed and modelled maximal radii (from some point within the profile) on the 
projection plane x = (£, rf) £ M 2 at direction angles aij (starting from a chosen coordinate direction for 
positive £, i] — 0): 

X^Et'max'KO-r^'K)] 2 - (14) 

ij 



() 



As the accuracy of is proportional to the size of the image, the sum (|14[) automatically takes this 

weighting into account (of course, direct weighting due to varying noise levels can be used as well) . 

We now represent the body B as a polytope [IT] . Let two vertices a and 6 of a facet have projection 
points (£ a ,Va)j {CbiVb), and (£,OiVo) be the point on the projection plane from which the radii and a are 
measured (this defines the profile offset that must be solved for in the inverse problem). With 

A = — sin a, B = cos a, 

C = r] a -r) h , L> = £ fc -£a, (15) 
E = A£ i0 + Br] Q , F = £ b ij a - £ Q ?/ h , 

the intersection point of the radius line and the projection of the facet edge ab is at 

DE- BF _ AF-CE 
K AD — BC V AD-BC { J 

and, to be in the correct direction of a and between the points a and &, the intersection point must satisfy 

(e-60(6-O>0, (£-a>)cosa>0, (17) 
{V -Va)(Vb — V) > (V - Vo) sin a > 0. 

If AD — BC = 0, the line in a-direction is parallel to the line ab, so there is no intersection unless the 
lines coincide, i.e., either of the numerators in (j!6[) vanishes. 

The model r max (a) can now be determined by going through all eligible facet edges and their intersection 
points x a b(a): 

^ax d) («) = max{||x o6 (a) - x \\\a,b e V+}, (18) 

where V+ is the set of vertices of the facets A+ that are both visible and illuminated. The set A+ of 
© is determined by ray-tracing [T3]. It is an approximation (correct to the order of the average facet 
area) of the actual visible and illuminated region, i.e., each facet either is or is not in A.+ (judging by 
its centroid): projection lines of obstructing facets inside a facet are neglected when the facets are small 
enough. 

The principle of using outer contours applies to AO data that do not generally show non-starlike or 
multiple contours (due to crater shadows) as the solar phase angles arccos(w, wq) are low and the reso- 
lution/deconvolution accuracy is not high. At high phase angles, even starlike bodies form non-starlike 
contours, and contours inside the outer contour appear in high-resolution images from probe flybys. 

The outer contour dO can be automatically derived for non-starlike shape models as well; such models 
can be constructed by, e.g., joining starlike submodels together, using a cylindrical coordinate frame 
or by determining the coordinates of a set of points with which a suitable surface (a new tessellation for 
each iteration) is defined via, e.g., mesh-free methods such as weighted/moving least squares |18) . For 
clarity, let us first assume that no other generalized profile contours exist outside dO. Denoting the edges 
of the facets of A+ by £ + , dO is constructed by the following algorithm: 

1. Construct the set IFq C £+ of the edges of £ + that are shared by a facet in A+ and by a facet not in 
A+ but for which (v, w) > 0. 

2. Construct the set T C £+ of the edges shared by a facet of A+ and a facet for which (i/, u) < 0. 

3. Construct the connected and ordered line sequences (lists of vertices) Eoi of the adjacent edges of J-$. 
The edges are defined by two vertices, and within Soi one vertex shares two edges. 

4. Construct the connected sequences from JF as in 3. 

5. The projections of the line sequences V(uj, Eoi) cannot intersect each other, but V(uj, Ej) can intersect 
each other and V(uj, Eoi) (intersection of projections can only occur when the surface folds away from 
sight, i.e., the line corresponds to a facet for which (u, lj) < 0). For any intersecting projected sequences, 
find the intersection points on the projection plane with the intersection test above. 
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6. Define the visible projections Eoi and Ej as the projected sequences V(u>, Eoi) and P(uj, E^) that may 
have a pij as an end point. 

7. Connect all Eoi and Ei that can form closed circuits (by systematically comparing the endpoints of 
the sequences). The circuit enclosing all the others (e.g., those due to shadows) is the approximation of 
the outer profile contour dO. 

If there are more than one generalized profile contours, the identification of the circuits should be arranged 
suitably to enable the comparison between the model and the data. For example, one shadow region 
inside dO and a smaller dOi outside dO due to a separate closed surface can be identified directly, and 
all circuits can be used in determining the best model. 

The position of a point in dO can be parametrized by using the path length along dO. The \% ^ s now 
given by (assuming one contour per profile) 

xl = IM^j) - Om(c i0 + cy) - x l0 )\\ 2 + A^(S* 4 - C 4 ) 2 , (19) 

ij i 

where o an m stand for observed and modelled, < c < 1 is the normalized path length along the 
measured and modelled contours dAii, dOi, Xio is the profile offset for each profile i, is the offset 
parameter for the path's starting point, A is a suitable weight factor, and 



Si = f ds, C, = f ds. (20) 
Thus, c = s/Si or c = s/Ci, where s is the usual path length. 

The contour fit can also be modelled by considering the distances of observed points from the model 
contour dO. Now we define 



xl 



= ^2mf{\\dO i (s) -xyll 2 }, (21) 



i j 



where Xij are the data points, and we label the points on dOi by s, and assume the translation due to 
Kio to be included in dOi(s); here it suffices to consider points in dO on whose normal lines xy lies. 
When dO is a set of line segments, the shortest distance required in (f2Tj) (let us denote it by <5 m j n ) is 
defined by: 

1. Let p'j £ R 2 be the projection of x on the line coinciding with the jth line segment (corresponding to 
tana = (£ a — Cb)/( 7 ?a — Vb) and (£0^0) ^ ( x -y) i n the intersection test above): 

-P 2 £ - CD? ] + CF , _ C 2 V - CDj + DF 

If the projection is inside the segment, let d'^ be the distance between and h. 

2. Let dk be the distance of x from the fcth end point of the line segments of dO. 
3- ^min is the smallest one of all the distances d'^ and dk. 

In addition to solving for the shape, we usually need to determine the target's spin state as well in order 
to have correct projection directions |12j . In most cases, the target revolves around a constant pole 
direction (f3, A) € S 2 at a constant rotation speed. The profile plane coordinates (£, -q) are the x' 2 - and 
Xg-components of 

x' = Rx, (23) 

where 

R = R y (^ - |)R*W> - A)R„(-j8)R,(-0o - n(* - to)), (24) 

where t is the time, f2 is the rotation speed (2ir/P for a constant rotation period P), <f>Q and the epoch 
to are some initial values, and ($,ip) € S 2 is the direction from the target to the observer. We determine 
(/3, A) and Q when solving the inverse problem. It is easy to accommodate other spin models such as 
precession |13j or nonconstant rotation speed \J3. in this formalism. 
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3.1 Regularization 



The parameters P describing the target usually have to be (moderately) regularized to prevent unrealistic 
solutions. One aspect is the smoothness of the body; the larger the target is, the less irregular it is expected 
to be. For some parts of the surface this is explicitly enforced by the profile data, so the regularization 
mostly pertains to the parts covered only by lightcurves that contain little information on nonconvex 
features. In those regions, undulation of the surface should be suppressed (the optimal choice of the 
suppression weight is discussed in section 4). 

For starlike bodies B, a simple (computationally x 2 -like) measure of global regularity is 

9S= f[r-(r)} 2 da, (25) 



JB 

where r is the radius of the model at the point corresponding to the surface element da; for polytopes, we 
can simply use gs = J2i( r i ~ ( r )) 2 - F° r such bodies, gs is typically quite efficient when the radius is given 
by a truncated Laplace series (itself a smoothing agent; see section 5 and [IT]) and the regularization 
weight is low. This is usually the case here as the profile contours already prevent runaway solutions, so 
gs is only needed to polish up the resolution level of shape detail. For higher weights or models with 
independent (uncorrelated) surface points, gs is not suitable as it emphasizes global roundedness more 
than local smoothness. 

A measure concentrating on local smoothness (and more suitable for more complex cases) can be con- 
structed by considering the negative values of the curvature function. For polytopes, a practical discrete 
version of this is computed by measuring how efficiently the facets not in the convex hull of the polytope 
can be blocked (from viewing or illumination) by their adjacent facets Taking into account the 

size and relative tilt angle of the possible blocker facets adjacent to the facet i, we can define, e.g., the 
following measure C by summing over the polytope and weighting suitably: 

^ = ^^^(1-^,^)), (26) 

where A% denotes the area of the facet i, and Aij the areas of those facets around it that are tilted above 
its plane |IIj (for i in the convex hull, Aij = by definition). In regularization, we minimize C (for 
convex bodies C = 0). 

A further smoothing constraint, to be used for non-star like contours dO if the observations do not cover 
the profile densely, is given by augmenting (1211) by 

A Et^/ m£ hdOi(8)-xf}ds, (27) 

which suppresses irregularity on surface parts not projected near the observed profile points. 

A physical constraint for most asteroids is that they are principal- axis rotators: their maximum moment 
of inertia is aligned with the rotation axis due to energy dissipation caused by the nonzero elasticity of 
the material of the body [30] . The regularization is defined by the symmetric inertia tensor [5] 




(28) 



where the inertia products P^ are 



Pij = / p(x)xiXj d 3 x, (29) 

JB 

and here we choose constant density p(x) = 1. We want to minimize the angle r between the z-axis of 
the model and the eigenvector / el 3 (normalized (/, I) — 1) corresponding to the largest eigenvalue of 
the inertia matrix I of the model shape B, so we can choose, for example: 

. g/ = (I-cos 2 T) 2 = [I-/ 3 (£) 2 ] 2 , (30) 
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where the square form if is useful for weighting purposes and for removing the sign ambiguity in 1$. A 
fast way of evaluating ^3 in (|30p for any polyhedron is described in [5] . Again, profile data constrain the 
result so strongly that usually the weight for gi is very low and sometimes can be set to zero to obtain, 
say, r < 4°. Enforcing a r much lower than this is seldom meaningful due to shape resolution level and 
inhomogeneities in the density. 

4 Optimal combination of data modes: maximum compatibility 
estimate 

From the statistical viewpoint, when we have two or more data modes, we consider their simultaneous 
probability distribution of model parameters and observations (augmented by prior or regularization 
distributions) in determining the posteriori distribution of the model and the corresponding estimates. 
The essential problem in this combining is inevitably the weighting of distributions. While the data 
modes share a common set of parameters describing the underlying model to be solved for, the models 
and modalities of observations may be completely different, and we seldom know a priori exactly how to 
compare and weigh their significance. 

Let us choose as goodness-of-fit measures (from which probability distributions can be constructed) the 
X 2 -functions of n data modes. Our task is to construct a joint xtot with well-defined weighting for each 
data mode: 

n 

xl t (P,D)=xl(P,D 1 )+J2^-iX-(P,D l ) D = {D i ,i = l,...,n} (31) 

i=2 

(to which regularization functions g(P) can be added), where Di denotes the data from the source i, and 
P € M. p is the set of model parameter values. We assume the x 2 -space to be nondegenerate, i.e., 

arg min \\ (P) ^ arg min (P), j 

In two dimensions, denote 

x(X) := { X ?|min Xt 2 ot ;A}, (32) 

vW '■= {x2l min xt 2 o t ; A l- 

The curve 

5(A) := [log x(A), log y(A)] (33) 

resembles the well-known "L-curve" related to, e.g., Tikhonov regularization PJ0I5]. However, here we 
make no assumptions on the shape of S. The curve S is a part of the boundary dlZ of the region 1Z £ R 2 
formed by the mapping x '■ K p — > R 2 from the parameter space P into x 2_ space: 

x = {p-> (log x 2 , log x')}, n = x {v) 

where the set V includes all the possible values of model parameters (assuming that x is continuous and 
well-behaved such that a connected 1Z and dlZ exist). If the possible values of xi are n °t bounded, the 
remaining part dlZ\S stretches droplet- like towards (00,00). The parameter A describes the position on 
the interesting part S C dlZ, and it is up to us to define a criterion for choosing the optimal value of A. 

The logarithm ensures that the shape of 5(A) is invariant under unit or scale transforms in the xf as they 
merely translate S in the (logx 2 ,logX2)"Pl anc - ^ a l so provides a meaningful metric for the logx 2 -space: 
distances depict the relative difference in x 2 -sense, removing the problem of comparing the absolute 
values of quite different types of xf- The endpoints of 5(A) are at A = and A = 00, i.e., at the values 
of x 2 that result from using only one of the data modes in inversion. We can translate the origin of the 
(logx?, logX2) _ pl ane to a more natural position by choosing the new coordinate axes to pass through 
these endpoints. Denote 

x = logx(A)| A=0 = logminx? (34) 
Vo = logy(A)| A ^.oo =logmmx2- 
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Then the "ideal point" (xo , yo) is the new origin in the (log x, log y)-plane. A natural choice for an optimal 
location on S is the point closest to (xo,yo), i.e., the parameter values Pq G P such that 



Pq = arg min 



([logxl(P) ~io? + [logx|(P) - y ] 2 ), 



(35) 



so we have, with A as argument, 



A = arg min ([log x(A) - x ] 2 + [logy (A) - y ? 



(36) 



In this approach, neither the numbers of data points in each xj nor the noise levels as such affect the 
solution for the optimal Pq as their scaling effects cancel out in each quadratic term. Pq is thus a pure 
compatibility estimate describing the best model compromise explaining the datasets of different modes 
simultaneously. 

We call the point Pq the maximum compatibility estimate (MCE), and Ao the maximum compatibility 
weight (MCW). This corresponds to the maximum likelihood estimate in the case of one data mode, 
or to the maximum a posteriori estimate as well since we can include regularization functions here. If 
regularizing is used, the weights for the functions are either determined in a similar manner (see below), 
or they can be fixed and the regularization terms are absorbed in \\ (otherwise S C dTZ does not hold). 

Another choice, frequently used in the L-curve approach, is to find the A at which S attains its maximum 
curvature [9l [7], but evaluating this point is less robust than finding Ao, and (|36|) is a more natural 
prescription, requiring no assumptions on the shape of S. We make two implicit assumptions here: 

1. The solutions Pon corresponding to points on dTZ should be continuous (and one-to-one) in P-space 
along dTZ at least in the vicinity of the solution corresponding to Ao- If this is not true (in practice, 
if Pa = arg min x^ ot (P) makes large jumps in P for various A around Ao), one should be cautious 
about the uniqueness and stability of the chosen solution P , and restrict the regions of P included 
in the analysis. 

2. The optimal point Ao on S should be feasible: if we have upper limits Ci to acceptable xh the feasible 
region T is the rectangle fliOogXi < logej. If [logxf (P ), logx|(P )] g T and T n TZ ^ 0, we 
choose the point on the portion 5cK closest to the one corresponding to Ao (i.e., logxf = l°g e ; 
for one i). If T n TZ = 0, the data modes do not allow a compatible joint model, so either the 
model is incorrect for one or both data modes, or one or both ei have been estimated too low (e.g., 
systematic errors have not been taken into account). Note that model insufficiency should be taken 
into account in the estimation of e,. 

Note that, in the interpretation TZ = x(P), \ Xtot an( ^ are au hi f &c t superfluous quantities, and 
we can locate the point estimate MCE Pq entirely without them with standard optimization procedures 
(and with no extra computational cost). However, it is useful (though computationally somewhat noisier) 
to approximate S via the minimization of Xtot with sample values of A (see Fig. 1), as in addition to 
obtaining the MCW Ao (and hence MCE as well) we can plot S to examine the mutual behaviour of the 
complementary data sources (including the position of the feasibility region J- w.r.t. S). The solution 
for Ao is also needed for constructing distributions based on Xtot- Another possibility to examine TZ and 
dJZ is direct adaptive Monte Carlo sampling, but this is computationally slow. 

This approach straightforwardly generalizes to n x 2 -functions and n — 1 parameters A^ describing the 
position on the n — 1-dimensional boundary surface dTZ of an ?i-dimcnsional domain TZ: the MCE is 




(37) 



and the MCW is 



A G W 



n — 1 



Aq = arg min 



i=l 



log 



Xj tot (A) l 2 
X?o ■ 



X 2 ,tot( A ) : = { 




(38) 
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Figure 1: S curve plotted for 2 Pallas with various weights A (LC for lightcurves, AO for adaptive optics 
profiles). 



Another scale invariant version of MCE can be constructed by plotting xf in units of xf /xfo and shifting 
the new origin to xf I xfo = 1 : 



Fo = argmm > g 1 ; A = argmm > g 1 ■ (39) 

L XiQ ■ - <■ Y.n 



2 4 . ^rx^A) 



This, however, is exactly the first-order approximation of (|3"T|) and (|55|) in 5 <C 1 when xf/xfo = 1 + 5, 
giving virtually the same result as (|37l) and (|38| as usually xf(Po)/xfo — 1 <IC 1 in the region around 
xf{Po), and any larger ratios of xf/xfo are n °t eligible for the optimal solution (see Fig. 1). 

Instead of the L2-norm x 2 (and the corresponding ^-distribution), we can choose some other goodness- 
of-fit measure e(P,D) > (and distribution) for the individual data modes. For a linear combination of 
these, we have 

n 

£ t ot(P, D) = £i(P, D{) + A 4 -l£ 4 (P, A)- 

»=2 

In lightcurve measurements, for example, the effect of systematic errors in both model and data dominates 
over random noise when the noise level is not high |15j . so it is not mandatory to use x 2 as a standard 
measure of fit. 

It is possible to use this approach for general regularizing functions g(P) as well (change xf — > g(P) for 
some i), but in such cases the shape of S must be taken into account. If it is possible to have a solution 
g(P') = for a regularizing function g (or an almost vanishing g(P') such that log g(P') — > — oo), the 
above scheme automatically returns P' and ignores the actual data altogether. Thus one should, e.g., set 
a lower practical limit to g(P) by looking at the shape of S, and choose the Ao within the restricted part 
of S. Likewise, one can use the above scheme for assigning noise- level- independent weights to subsets of 
the same data mode (rather than have the standard x 2 evaluated from all data points), but obviously 
the subsets cannot be chosen arbitrarily if the result is to make sense. For example, one can estimate the 
optimal weight for one lightcurve that appears to reveal features not contained in other lightcurves and 
thus judge its real significance. Even one noisy lightcurve with a few points, taken at a special observing 
geometry, may well contain significant information that needs to be weighed more against less noisy but 
more ordinary lightcurves. 
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Figure 2: Sample observed (solid lines) vs. modelled (dashed lines) AO contours for 2 Pallas. Coordinates 
are in pixel units. 



5 Numerical implementation 

As examples of the optimal combining of lightcurves and AO profiles, we show representative results for 
the asteroids 2 Pallas and 41 Daphne. Full detailed descriptions of the observations and models of these 
targets are presented in [3] and Carry et al. (in preparation). An example of an even more irregular shape 
constructed with our procedure is the model of the primary body of the binary asteroid 121 Hermione 0]. 
The lightcurve x\ was computed as in [IT] [12] and profile x% as m ^ ne starlike case of section 3, and the 
minimization of Xtot was performed as in |11[ 112] . The observed profiles are projections of the target on 
the plane-of-sky S 2 converted to pixels on the instrument plane, while the model is constructed in absolute 
(km) size, so the model/profile scale conversion is given by the AO instrument's angular resolution and 
the distance between the target and the observer. The profile contour extraction procedure with wavelets 
(as an average of several AO images obtained in a short time interval) is described in [2j |3] . 

In general, the resolution of the model must be somewhat lower than the apparent resolution of the AO 
images as the sparse profile samples will produce artificial features elsewhere in the model if a near-perfect 
profile fit is enforced (even if the observed profile details were exactly right). The inverse problem has 
thus some ill-posedness at local scales starting near the profile resolution level, but the ill-posedness at 
more global scales, inherent to lightcurve data [12] [14] [15] , is removed with AO profiles. The weight factor 
A mostly takes care of this, and fine-tuning is obtained with As for the smoothness constraint gg. For the 
examples here, the weight of the inertia regularization function gi was low as there were several profiles 
available; virtually the same result was achieved with A/ = 0. The weights A and As were determined 
with the scheme of section 4; the examined interval of As was restricted to realistic values corresponding 
to the resolution level of the AO images. 

Fig. 1 depicts a typical evaluation of the curve S at various choices of A; or rather, this plot portrays the 
cross-section of the 2-surface dlZ in M 3 with As fixed at its final optimal value. The values for xf are 
normalized to be the rms deviations of model fits di — y/ xf /Ni , as in logarithmic scale this corresponds 
only to a shift of origin and a uniform linear change of plot scaling. The plotted points outline the curve 
5(A) that is rather an oblique line than an L-shape, and the ideal point region, i.e., the point closest to the 
lower left-hand corner, can directly be found. The endpoints A = and A = oo stop at saturation regions 
rather than continue to large distances in the logx 2 -space. As can be seen from Fig. 1, computational 
noise in the estimated points at A = and A = oo, corresponding to a small change of the position of the 
new origin w.r.t. S, does not affect the estimated location of the optimal point on S significantly. 

Sample observed vs. modelled profiles for 2 Pallas and 41 Daphne are shown in Figs. 2 and 3. The starlike 
surface model was described by the exponential Laplace (spherical harmonics) series for the surface radius 



truncated at suitable I, rn, with c; m as the shape parameters to be solved for. Other model parameters are 



r m 




(40) 



I I II 
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Figure 3: Sample observed (solid lines) vs. modelled (dashed lines) AO contours for 41 Daphne. Coor- 
dinates are in pixel units. 



the profile offset (£o, Vo) f° r each image and the spin parameters. For asteroid 2 Pallas (a rather spherical 
body with size class 500 km), the Laplace series was truncated at maximal I = 6, m = 6, while for the 
more irregular 41 Daphne (size class 200 km) the truncation point I = 8, to = 6 was more appropriate. 
The early truncated Laplace series and the choice of the truncation point are implicit regularization 
measures as such. We leave the discussion of the choice of model discretization level elsewhere (cf. |17j ) 
as here its effect on the data mode weighting was negligible (within a feasible set of choices), and the 
resolution level of AO images (as well as keeping As low and avoiding artificial surface features) essentially 
determined the choice in practice after some sampling. For AO data, the choice of the Laplace series as 
a model is practical, while for, e.g., detailed space probe data a mesh of independent surface points is 
more accurate and computationally feasible. 

Once the weight factors A and As are determined, the result is usually stable and restricted to one region 
in the parameter space P: probing feasible solutions P corresponding to Xtot(P) slightly lower than 
Xtot(-Po) produces essentially the same results. Due to restricted orbital geometries, lightcurve data alone 
often imply two almost equally possible pole directions with mirror- like shape solutions [14] [15] ; even one 
AO (or other) image usually resolves this typical ambiguity [19]. The result is also typically stable w.r.t. 
weights in the vicinity of MCW. The obtained MCE appears to be well justified when one samples the 
solutions along S: it provides a very good match to profile details without straying far from the observed 
lightcurves, and does not predict too prominent features on the parts of the surface not projected onto 
the profile contours. 

6 Conclusions and discussion 

We have examined the classes of shapes reconstructable by the (generalized) profiles of objects in M 3 , and 
presented a method for using lightcurves and the observed contours of generalized profiles simultaneously 
to produce shape (and spin) models with more details (and a lower degree of ill-posedness) than in the 
pure lightcurve mode. We have also shown that there is a well-justified criterion and an efficient method 
for determining the optimal weighting of data modes. Applied to real data, the method works very well, 
and we can use simple regularization functions. In addition to adaptive optics observations, asteroid 
profiles can also be obtained from other sources such as interferometry, space telescopes, and stellar 
occultations (partial profiles). 

The use of profiles is practical as it removes two sources of systematic errors inherent to using full images 
(brightness distributions I on the image plane) : the errors in I from AO deconvolution and the model 
I errors due to the insufficently modellable light-scattering properties of the surface of the target body. 
On the other hand, profile determination requires the data to be sharp enough, not with fuzzy images. 
If the images are fuzzy, we usually have to resort to using some brightness and blurring model for fitting 
full images, even though the result will be less certain. 
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The concept of the maximum compatibility estimate is directly applicable to any inverse problems with 
complementary data modes. The invariance properties of the MCE make it more generally usable than 
heuristic strategies for choosing the weights, especially when they use assumptions on the shape of dlZ 
or other case-specific characteristics. 

Acknowledgements 

It is a pleasure to thank Benoit Carry and Josef Durech for discussions and comments. The sample 
adaptive optics data used in figures here are courtesy of B. Carry, A. Conrad, J. Drummond, C. Dumas, S. 
Erard, and W. Merline. This work was supported by the Academy of Finland (project "New mathematical 
methods in planetary and galactic research" ) . 

References 

[1] M. Beige, M. Kilmer, and E. Miller, Efficient determination of multiple regularization parameters in 
a generalized L-curve framework, Inverse Problems, 18 (2002), 1161-1183. 

[2] B. Carry, C. Dumas, M. Fulchignoni, W. Merline, J. Berthier, D. Hestroffcr, T. Fusco, and P. Tam- 
blyn, Near-infrared mapping and physical properties of the dwarf-planet Ceres, Astron. Astrophys., 
478 (2008), 235-244. 

[3] B. Carry, C. Dumas, M. Kaasalainen, and 9 colleagues, Physical properties of 2 Pallas, Icarus, (2009) 
in press. 

[4] P. Descamps and 22 colleagues, New insights on the binary asteroid 121 Hermione, Icarus, 203 
(2009), 88-101. 

[5] A. Dobrovolskis, Inertia of any polyhedron, Icarus, 124 (1996), 698-704. 

[6] J. Durech and M. Kaasalainen. Photometric signatures of highly nonconvex and binary asteroids, 
Astron. Astrophys., 404 (2003), 709-714. 

[7] H. Engl and W. Grever, Using the L-curve for determining optimal regularization parameters, Numer. 
Math., 69 (1994), 25-31. 

[8] H. Goldstein, "Classical mechanics" (second edition), Addison- Wesley, Reading, Mass., 1980. 

[9] M. Hanke, Limitations of the L-curve method in ill-posed problems, BIT, 36 (1996), 287-301. 

[10] M. Kaasalainen, L. Lamberg, K. Lumme, and E. Bowell, Interpretation of lightcurves of atmosphere- 
less bodies. I. General theory and new inversion schemes, Astron. Astrophys., 259 (1992), 318-332. 

[11] M. Kaasalainen and J. Torppa, Optimization methods for asteroid lightcurve inversion. I. Shape 
determination, Icarus, 153 (2001), 24-36. 

[12] M. Kaasalainen, J. Torppa, and K. Muinoncn, Optimization methods for asteroid lightcurve inver- 
sion. II. The complete inverse problem, Icarus, 153 (2001), 37-51. 

[13] M. Kaasalainen, Interpretation of lightcurves of precessing asteroids, Astron. Astrophys., 376 (2001), 
302-309. 

[14] M. Kaasalainen and L. Lamberg, Inverse problems of generalized projection operators, Inverse Prob- 
lems 22 (2006), 749-769. 

[15] M. Kaasalainen and J. Durech, Inverse problems of NEO photometry: Imaging the NEO population, 
in "Proceedings of IAU: Symposium 236", 2, Milani, Valsecchi, and Vokrouhlicky, eds., Cambridge 
(2007), 151-166. 



15 



[16] M. Kaasalainen, J. Durech, B. Warner, Y. Krugly, and N. Gaftonyuk, Acceleration of the rotation 
of asteroid 1862 Apollo by radiation torques, Nature, 446 (2007), 420-422. 



[17] J. Kaipio and E. Somersalo, "Statistical and computational inverse problems", Springer, New York 
2005. 

[18] D. Levin, The approximation power of moving least squares. Math. Comp., 67 (1998), 1517-1531. 

[19] F. Marchis, M. Kaasalainen, E. Horn, J. Berthier, J. Enriquez, D. Hestroffer, D. Lc Mignant, and I. 
de Pater, Shape, size and multiplicity of main-belt asteroids. I. Keck adaptive optics survey, Icarus, 
185 (2006), 39-63. 

[20] P. Pravec, A. Harris, and T. Michalowski. Asteroid rotations, in "Asteroids III", Bottke, Cellino, 
Paolicchi, and Binzel, eds., U. Arizona Press, Tucson (2002), 113-122. 

[21] S. Savarese, M. Andretto, H. Rushmeier, F. Bernardini, and P. Perona, 3D Reconstruction by Shadow 
Carving: Theory and Practical Evaluation, International Journal of Computer Vision, 71 (2007), 
305-336. 

E-mail address: First. Lastname [at] tut.fi 



16 



